/* 
 * File:   main.cpp
 * Author: Cristiano Strieder // cstrieder@gmail.com
 *
 * Created on April 22, 2009, 11:32 AM
 */

#include <stdlib.h>
#include <iostream>
#include <string>
#include <fstream>
#include <sstream>
#include <vector>
#include <math.h>

using namespace std;

int lines, cols; // Store the size(lines and cols) of the matrix

/*
 * 
 */

int matrix_lines(istream &iFile) {
// Método para obter o número linhas da matriz.
    iFile.clear();
    iFile.seekg(0);
    int lines = 0;
    string line;
    while(!iFile.eof()) {
        getline (iFile, line);
        lines++;
    }
//    cout << "Número de linhas: " << lines-1 << endl;
    return lines-1;
}

int matrix_cols(istream &iFile) {
// Método para obter o número de colunas da matriz.
    iFile.clear();
    iFile.seekg(0);

    string line;
    double n;
    int cols = 0;

    getline (iFile, line);
//    cout << line << endl;

    stringstream is(line);
    n = 0;
    while (is >> n) {
//        cout << n << endl;
        cols++;
    }


//
//    cout << "Número de colunas: " << cols << endl;
    return cols;
}

vector<vector<double> > matrixRead(istream &iFile) {
// Método para ler a matriz.

    vector<vector<double> > matrix( cols , vector<double>( lines ) );

    double value;
    iFile.clear();
    iFile.seekg(0);
    
    for(int i = 0; i < lines; i++){
        for(int j = 0; j < cols; j++){
           matrix[i][j]  = (iFile >> value,value);
           cout << matrix[i][j] << " ";
        }
        cout << "\n";
    }

    return matrix;
}

vector<vector<double> > gradX(vector<vector<double> > matrix) {

    vector<vector<double> > grad_x( cols , vector<double>( lines ) );

    for(int i = 0; i < lines; i++) {
        for(int j = 0; j < cols; j++) {

            if(j==0) { // primeira coluna
                grad_x[i][j] = matrix[i][1] - matrix[i][0];
            }else if(j==cols-1) { // ultima coluna
                grad_x[i][j] = matrix[i][j] - matrix[i][j-1];
//                cout << "matrix[i][j-1]= " << matrix[i][j-1] << " matrix[i][j]= " << matrix[i][j] << endl;
            }else // restante da matriz
                grad_x[i][j] = (matrix[i][j+1] - matrix[i][j-1])/2;
        }
    }
    cout << "gradiente X" << endl;

    for(int i = 0; i < lines; i++){  // Output gradient X
        for(int j = 0; j < cols; j++){
           cout << grad_x[i][j] << " ";
        }
        cout << "\n";
    }

    return grad_x;
}

vector<vector<double> > gradY(vector<vector<double> > matrix) {

    vector<vector<double> > grad_y( cols , vector<double>( lines ) );

    for(int i = 0; i < lines; i++) {
        for(int j = 0; j < cols; j++) {

            if(i==0) { // primeira linha
                grad_y[i][j] = matrix[1][j] - matrix[0][j];
            }else if(i==lines-1) { // ultima linha
                grad_y[i][j] = matrix[i][j] - matrix[i-1][j];
            }else // restante da matriz
                grad_y[i][j] = (matrix[i+1][j] - matrix[i-1][j])/2;
        }
    }
    cout << "gradiente Y" << endl;

    for(int i = 0; i < lines; i++){  // Output gradient Y
        for(int j = 0; j < cols; j++){
           cout << grad_y[i][j] << " ";
        }
        cout << "\n";
    }

    return grad_y;
}

vector<vector<double> > matrixNorm(vector<vector<double> > matrix) {

    vector<vector<double> > matrix_norm( cols , vector<double>( lines ) );
    double max = 0;
    int max_col,max_line;

    for(int i = 0; i < lines; i++){  // Encontra o máximo da matrix
        for(int j = 0; j < cols; j++){

           if(matrix[i][j] > max){
               max = matrix[i][j];
               max_col = j;
               max_line = i;
           }

        }
    }
    cout << "Max value " << max << " at line " << max_line+1 << ", col " << max_col+1 <<  endl;

    for(int i = 0; i < lines; i++){  // Normaliza matrix
        for(int j = 0; j < cols; j++){

            matrix_norm[i][j] = matrix[i][j]/max;

        }
    }
    cout << "Matrix normalizada" << endl;

    for(int i = 0; i < lines; i++){  // Output gradient Y
        for(int j = 0; j < cols; j++){
           cout << matrix_norm[i][j] << " ";
        }
        cout << "\n";
    }

    return matrix_norm;
}

class vectors {
public:
    double module, phase;
};

vector<vector<double> > matrixModPhase(
                                        vector<vector<double> > grad_x,
                                        vector<vector<double> > grad_y
                                      ) {

    vectors matrix_mod_phase[lines][cols]; // matrix de modulos e phases
    double module,phase;

    for(int i = 0; i < lines; i++){
        for(int j = 0; j < cols; j++){
            
            module = sqrt((double)pow(grad_x[i][j],2)+(double)pow(grad_y[i][j],2));
            cout << module << " ";

            phase = atan(grad_y[i][j]/grad_x[i][j]);
            cout << phase << endl;

            matrix_mod_phase[i][j].module = module;
            matrix_mod_phase[i][j].phase = phase;
        }
        cout << "\n";
    }

    vector<vector<vector<double> > > matrix_mod_phase1( cols , vector<double>( lines ), vector<vector<double>( 2 ) );
    // The last dimension store:
    //                      1 - module
    //                      2 - phase

    for(int i = 0; i < lines; i++){  // Output gradient Y
        for(int j = 0; j < cols; j++){
           cout << matrix_mod_phase[i][j].module << " ";
           cout << matrix_mod_phase[i][j].phase << " ";
        }
        cout << "\n";
    }

    vector<vector<double> > > matrix_mp;
    return matrix_mp;
}

int main(int argc, char** argv) {
    
    
    string filename;
    filename = "/home/cristiano/matrix4.dat";
    ifstream iFile;
    iFile.open(filename.c_str());

    if(!iFile) { // file couldn't be opened
        cerr << "Error: file could not be opened" << endl;
        exit(1);
    }

    cols = matrix_cols(iFile); // Get matrix size
    lines = matrix_lines(iFile);
    cout << "Matrix " << cols << "x" << lines << " - " << filename << endl;

    vector<vector<double> > matrix;

    matrix = matrixRead(iFile); // Load matrix

    vector<vector<double> > matrix_norm;

    matrix_norm = matrixNorm(matrix); // Normaliza matrix

    vector<vector<double> > grad_x,grad_y;

    grad_x = gradX(matrix_norm); // Calc gradient
    grad_y = gradY(matrix_norm); // Calc gradient

    vector<vector<double> > matrix_mod_phase;

    matrix_mod_phase = matrixModPhase(grad_x, grad_y); // Calc Modulo e fase

    iFile.close();
// Opa2
    return (EXIT_SUCCESS);
}

